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ABSTRACT 

Using cosmological A^-body simulations we study the radial velocity distribution 
in dark matter haloes focusing on the lowest-order even moments, dispersion and 
kurtosis. We determine the properties of ten massive haloes in the simulation box 
approximating their density distribution by the NFW formula characterized by the 
virial mass and concentration. We also calculate the velocity anisotropy parameter 
of the haloes and find it mildly radial and increasing with distance from the halo 
centre. The radial velocity dispersion of the haloes shows a characteristic profile with 
a maximum, while the radial kurtosis profile decreases with distance starting from a 
value close to Gaussian near the centre. We therefore confirm that dark matter haloes 
possess intrinsically non-Gaussian, flat-topped velocity distributions. We find that the 
radial velocity moments of the simulated haloes are quite well reproduced by the 
solutions of the Jeans equations obtained for the halo parameters with the anisotropy 
measured in the simulations. We also study the radial velocity moments for a composite 
cluster made of ten haloes out to ten virial radii. In this region the velocity dispersion 
decreases systematically to reach the value of the background, while kurtosis increases 
from below to above the Gaussian value of 3 signifying a transition from a flat-topped 
to a strongly peaked velocity distribution with respect to the Gaussian, which can be 
interpreted as the dominance of ordered flow with a small dispersion. We illustrate 
the transition by showing explicitly the velocity distribution of the composite cluster 
in a few radial bins. 

Key words: methods: A-body simulations - methods: analytical - galaxies: clusters: 
general - galaxies: kinematics and dynamics - cosmology: dark matter 


1 INTRODUCTION 

The density distribution of bonnd structnres has been a sub¬ 
ject of vigorous research in recent years. An NFW (Navarro, 
Frenk & White 1997) density profile was established as a uni¬ 
versal formula describing dark matter haloes in large mass 
range arising in cosmological A-body simulations (for re¬ 
cent refinements see Navarro et al. 2004; Diemand, Moore 
& Stadel 2004a; Tasitsiomi et al. 2004; Merritt et al. 2005). 
Much less attention has been devoted to the velocity distri¬ 
bution in bound structures, although it also carries substan¬ 
tial amount of information about the formation and evolu¬ 
tion of structure in the Universe. The moments of the ve¬ 
locity distribution in particular proved extremely useful in 
studies of mass distribution in galaxies and clusters (Kro- 
nawitter et al. 2000; van der Marel et al. 2000; Lokas & 
Mamon 2003). 


Recently a few groups studied the velocity distribution 
of simulated haloes showing, for example, that it departs 
significantly from a Gaussian distribution expected for viri- 
alized objects (Kazantzidis et al. 2004; Sanchis, Lokas & 
Mamon 2004; Diemand Moore & Stadel 2004b). Such be¬ 
haviour, manifesting itself in a form of kurtosis value lower 
than 3, was also seen in real objects like the Coma cluster 
(Lokas & Mamon 2003) or the Draco dwarf (Lokas, Mamon 
& Prada 2004). There have been few attempts to provide a 
theoretical explanation for departures from Gaussianity; it 
turns out however that weakly non-Gaussian velocity distri¬ 
butions arise naturally in Jeans theory of equilibrium struc¬ 
tures (Merrifield & Kent 1990; Lokas 2002; Kazantzidis et 
al. 2004). Hansen et al. (2004) have shown that the flat¬ 
tened velocity distribution can also be interpreted in terms 
of Tsallis statistics. 
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In this Letter we study the radial velocity distribution 
of dark matter haloes resulting from cosmological A^-body 
simulations focusing on the most useful even moments, the 
dispersion and kurtosis. We show that they are quite well 
reproduced by the solutions of the Jeans equations. We also 
look at the behaviour of the moments outside the virial ra¬ 
dius and the overall velocity distribution. 

2 THE SIMULATED DARK MATTER HALOES 

For this study we run a cosmological dark matter sim¬ 
ulation within a box of size 150 h~^ Mpc assuming the 
concordance cosmological model (ACDM) with parameters 
IJm = 0.3, Ua = 0.7, h = 0.7 and as = 0.9. We have used 
a new MPI (Message-Passing Interface) version of the origi¬ 
nal ART (Adaptive Refinement Tree) code (Kravtsov et al. 
1997). The ART code achieves high spatial resolution by re¬ 
fining the underlying uniform grid in all high-density regions 
with an automated refinement/derefinement algorithm. In 
the new MPI version of the code the whole simulation box 
is subdivided into cuboids which contain approximately the 
same numerical tasks to achieve load balance (not necessar¬ 
ily the same number of particles). Each cuboid is handled 
by one node of the computer using internal OpenMP par¬ 
allelization. At the given node the remaining box is evolved 
using more massive particles as described for the multi-mass 
version of the ART code (Klypin et al. 2001). This procedure 
minimizes the communication between nodes. After each ba¬ 
sic integration step the borders of the cuboids can be moved 
if load balance changes. The simulation used 256® particles, 
therefore it achieved a mass resolution of 1.7 x 10®° 

The basic grid of the simulation was 256, the maximum re- 
hnement level was 6. Therefore, the force resolution (2 cells) 
reached 18/i“® kpc. 

The halos have been identified using the hierarchical 
friends-of-friends (HFOF) algorithm (Klypin et al. 1999) us¬ 
ing different linking lengths. The linking length of 6 = 0.17 
extracts haloes with masses close to the virial mass for our 
cosmological model. (We define the virial mass and radius as 
those with mean density Ac = 101.9 times the critical den¬ 
sity, see Lokas & Hoffman 2001.) However, the FOF method 
is known to connect particles along thin bridges which are 
likely to be broken once the linking parameter is changed. 
This also means that the shape of structure formed by linked 
particles differs strongly from spherical and the centre of a 
halo does not necessarily coincide with the region of maxi¬ 
mum density. 

In order to find points of maximum density we cut out 
the particles out to a virial radius (assuming that the FOF 
mass is equal to the virial mass) and calculate the centres 
of mass in spheres of decreasing radius, each time making 
the centre of mass the new halo centre. The centres deter¬ 
mined in this way agree very well with centres found with 
much smaller linking lengths which pick up regions of higher 
density. We then determine the new masses and virial radii 
of the haloes by summing up masses of all particles found 
inside the virial radius defined as before but now measured 
with respect to the new centre. Choosing the haloes for the 
analysis we have generally taken the most massive objects 
found in the simulation box. We have, however, rejected 
those undergoing a major merger with two or more sub¬ 
haloes of comparable mass which would significantly depart 


Table 1. Properties of the simulated haloes. 


Halo 


Tv 

[Mpc] 

c 

1 

7.49 

2.33 

5.968 

2 

7.04 

2.29 

10.098 

3 

6.83 

2.26 

5.360 

4 

6.49 

2.23 

8.301 

5 

5.87 

2.15 

6.285 

6 

5.35 

2.08 

10.320 

7 

4.82 

2.02 

8.577 

8 

4.25 

1.93 

3.552 

9 

3.86 

1.87 

7.607 

10 

3.54 

1.82 

7.424 


from equilibrium. Within the virial radius, the haloes typi¬ 
cally have a few times 10® particles. The virial masses and 
radii of the haloes, and r„, are listed in Tabled 

For each of the ten haloes we determined its den¬ 
sity profile taking averages in radial bins of equal loga¬ 
rithmic length up to the virial radius. The data points ob¬ 
tained in this way were assigned errors estimated as Pois¬ 
son fluctuations. The measured density profiles are well 
approximated by the NFW formula (Navarro et al. 1997) 
e(s)/ec,o = Ac c ^ 5 ( c )/[3 s (1-I-cs)^] where s = rjr^, £ic,o is 
the present critical density, Ac = 101.9 is the characteris¬ 
tic density parameter, c is the concentration parameter and 
g{c) = [ln(l -I- c) — c/(l -I- c)]“®. Fitting this formula to our 
measured density profiles we determined the concentration 
parameters of the haloes which are listed in Tabled Our es¬ 
timated concentrations are consistent with the dependence 
of c on mass inferred from A-body simulations by Bullock 
et al. (2001), also run with a ACDM cosmology. 

In the following, all velocities will be calculated with re¬ 
spect to the mean velocity inside the virial radius of a given 
halo. For each halo we measure the mean radial velocity 
(with the Hubble flow added) of dark matter particles en¬ 
closed in shells of thickness 0.1r„ centered at 0.05r„, 0.15r„ 
etc. assuming a convention that the negative sign indicates 
infall motion towards the centre of the halo and normaliz¬ 
ing to Vv = {GMv, the circular velocity at r„. For 
a fully virialized object the mean radial velocity should be 
zero. We find that the measured values are consistent on 
average with zero inside for most haloes, although some 
radial variations caused by internal streaming motions are 
present. They do not exceed 0.4Vj, however, indicating a 
rather relaxed state of the objects. 

We then average our measurements calculating the 
mean and dispersion of the measurements for the ten haloes 
in the corresponding bins. The results are plotted in the up¬ 
per left panel of Fig. Q with thicker solid line showing the 
mean and the two thinner lines indicating the dispersion 
with respect to the mean. The upper right panel of Fig. 0 
shows the radial profile of the anisotropy parameter 


averaged in the same way over the ten haloes, where af = 
ag + and ag^^ and Cr are the angular and radial velocity 
dispersions with respect to the mean velocities. As we can 
see, the orbits are mildly radial, with positive mean /3. 
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Figure 1. The mean radial velocity in units of Vv, the circular ve¬ 
locity at Vv (upper left panel), the radial profile of the anisotropy 
parameter [3 (upper right panel) and the radial moments: veloc¬ 
ity dispersion ar in units of Vv (lower left panel) and kurtosis 
Kr (lower right panel) inside Vv averaged over ten haloes listed 
in Table Q Thicker solid lines show the mean and the thinner 
solid lines the dispersion in the measurements. Dashed lines give 
predictions from the Jeans equations for c = 7.4 and a variable 0 
changing as shown in the upper right panel. 


The dispersions needed for the calculation of /3 were ap¬ 
proximated hy a = S where is the most natural estimator 
of the variance from a large sample of n velocities Vi 



i=l 


where v = Vi)/n is the mean of velocities in the shell. 

The lower left panel of Fig. Q shows the radial velocity dis¬ 
persion calculated in the same way. In the lower right panel 
of Fig.0we plot the fourth radial velocity moment normal¬ 
ized by the variance squared, i.e. the kurtosis 


K = 


(5'2)2 


(3) 


Although the moments show a significant scatter from 
halo to halo and some variability due to substructure in 
every halo, the overall trend in their behaviour, shown by 
the mean values plotted in Fig. Q is clearly visible. In par¬ 
ticular, the slightly increasing and then decreasing velocity 
dispersion profile agrees well with the predictions of Lokas 
& Mamon (2001) for the NFW haloes. The trend observed 
in the behaviour of kurtosis which decreases from a value 
close to 3 (characteristic of a Gaussian distribution) to val¬ 
ues well below 3 suggests a flat-topped velocity distribution 
in the outer radial bins. 


3 COMPARISON WITH PREDICTIONS FROM 
THE JEANS FORMALISM 

The Jeans formalism (e.g. Binney & Tremaine 1987) relates 
the velocity moments of a gravitationally bound object to 
the underlying mass distribution. We summarize here the 
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formalism, as developed in Lokas (2002) and Lokas & Ma¬ 
mon (2003). The second and fourth-order vf. radial ve¬ 
locity moments obey the Jeans equations 


d / 


2^ , 2/3 2 

) 3- 


_d 

dr 


■^{vVr) + —uvf + Svcr^^f- 


d$ 

dr 


0 

0 


(4) 

(5) 


where i/ is the 3D density distribution of the tracer popu¬ 
lation (here it is the same as the total mass density) and 
<1? is the gravitational potential, which for an NFW density 
distribution is ^{s)/Vf = —g(c)ln(l -|- cs)/s. The second 
equation was derived assuming the distribution function of 
the form f{E,L) = fo{E)L~^^ and the anisotropy param¬ 
eter /3, equation 0 , to be constant with radius. We will 
consider here — oo < /3 < 1 which covers all interesting pos¬ 
sibilities from radial orbits (/3 = 1) to isotropy (/3 = 0) and 
circular orbits (/3 —> —oo). 

The solutions to equations for /3 = const are 


u{r)a^{r) = r f x^^zz(a;)^5- da; (6) 

J r 

v{r)vRr) = dx . (7) 

Jr da: 

After introducing expression into 0 and inverting the 
order of integration in 0 the expression for the fourth mo¬ 
ment reduces to a single integral 


u{r)vf{r) = 3r 


-20 


..2/3 


d$ 

v{x)^—[^{x) — ‘I>(r)] da: . (8) 

da; 


For the specific case of NFW-distributed dark matter 
particles tracing their own gravitational potential the formu¬ 
lae reduce to the following expressions which can be easily 
calculated numerically 


l{s) 


VJ 


vRs) 

VJ 


^^{l + csfg{c) 

^20-3 


I 


ln(l + cz) — 


1 -I- cz. 


(1 -I- C2:)2 ; 
3s^“^^(l + cs)^g^{c) 

ln(l-h cs) ln(l-hC2:) 


dz 


(9) 


( 10 ) 


1 -I- cz. 


dz 


where the distances were expressed in units of the virial 
radius (s = r/r„, 2 = a;/r„) and Vv is the circular velocity 
at the virial radius. As usual, we express the fourth-order 
moment in terms of the dimensionless radial kurtosis 


Kr{r) 


Vrjr) 

at{r)' 


( 11 ) 


Although the velocity dispersion profile depends on the 
concentration parameter (see Lokas & Mamon 2001), our 
simulated haloes all have similar concentration, so for the 
comparison with the predictions from the Jeans formalism 
we will adopt their mean concentration of c = 7.4 and focus 
on the dependence on the anisotropy parameter /3. Since the 
anisotropy changes roughly between /3 = 0 in the centre and 
/3 = 0.6 at the virial radius, as shown in the upper right 
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r/r^ r/r^ 

Figure 2. Predicted radial velocity moments for c = 7.4 and 
different 0. The left (right) panel shows the radial velocity dis¬ 
persion (7r (kurtosis Kr) for 0 = 0 with the lower (upper) solid 
line and for 0 = 0.6 with the upper (lower) solid line. The dashed 
line in each panel shows the result for 0 linearly increasing from 
0 to 0.6. 

panel of Fig. Q we first plot in Fig. 121 the velocity moments 
obtained from the solutions a-iTnii for these two constant 
values of f}. We do not have to restrict our analysis to a 
constant /3 however, since the main contribution to ar{r) 
and Kr (r) comes from the region close to r in the integrands 
given in 0-0 and the moments are actually sensitive only 
to a local value of /3 at a given r. We can therefore reproduce 
the moments for variable beta taking a different (constant) 
j3 for each r. Assuming that f3 grows linearly between 0 
and 0.6 we get results shown with a dashed line in each of 
the panels in Fig. |21 We see that the new solution changes 
smoothly from the result for /3 = 0 near the centre to the 
result for jd = 0.6 near the virial radius. 

In the lower panels of Fig. 0 we show with dashed lines 
the predictions from R-fTil for (3 changing as shown in the 
upper right panel. Namely, as described above, for the cal¬ 
culation of the predicted moments in each bin we assume a 
value of /3 found in the same bin after averaging over ten 
haloes. We can see that the solutions of the Jeans equations 
reproduce quite well the trends in the behaviour of the ra¬ 
dial velocity moments, especially in the case of dispersion. 
Although the solutions are not formally exact, they do pro¬ 
vide fairly good approximations. 

4 RADIAL MOMENTS BEYOND THE VIRIAL 

RADIUS 

In this section we extend the analysis beyond the virial ra¬ 
dius and look at the behaviour of the moments there. Fig.|21 
again shows the same quantities as in Fig. Q but for a com¬ 
posite cluster made of particles from the ten haloes with 
distances scaled to of each halo. We plot the measure¬ 
ments up to distances of 10 r„, averaged in radial bins of 
size 0.2r„. The mean radial velocity displays a characteristic 
behaviour: just outside the virial radius it drops to negative 
velocities signifying infall, reaches the turn-around radius 
of the cluster at about 3-4 r„ and approaches the Hubble 
flow at large distances. The velocity dispersion falls to the 
value characteristic of the background just outside 2 r„. The 
kurtosis falls down inside r„ to values below 3 and then in¬ 
creases sharply at larger distances and oscillates around 5 
but staying well above 3. This signifies a velocity distribution 
which is more peaked in the centre than the Gaussian. The 
solid line in the upper right panel of the Figure shows the 
anisotropy parameter defined as in eq. 0 . For comparison 
we also plot with a dashed line the parameter with the con- 



Figure 3. The same quantities as in Fig. ^measured for a com¬ 
posite cluster made of ten haloes for distances out to 10 r„. The 
dashed line in the upper right panel shows the 0 parameter in¬ 
cluding the non-zero mean velocities (see text). 


tribution from non-zero mean velocities (i.e. R = 1 —Uj/2ur 
where vf = + vi^). The anisotropies drop to values below 

zero outside r„ and approach zero or unity at large distances. 
This behaviour can be understood as due to dominance of 
ordered flow in the radial direction outside the virial radius: 
random motions in the radial direction are smaller while the 
tangential ones remain the same which results in negative j3. 


5 SUMMARY AND CONCLUSIONS 

We summarize our results by plotting explicitly the radial 
velocity distribution of our composite cluster in Fig.|2] The 
histograms shown in the Figure were made for radial bins of 
equal number (4 x 10"*) of particles and show their velocity 
distribution with respect to the mean velocity in each bin. 
The range of particle positions in each bin (increasing from 
top to bottom panel) is marked in the corresponding panel. 
On top of the histograms we plot with dashed lines the nor¬ 
malized Gaussian distributions with dispersions as shown in 
the lower left panel of Fig.|21 

The first two panels show the distribution inside the 
virial radius. We can see that the distribution is Gaus¬ 
sian to a very good approximation only in the very cen¬ 
tre. In the next bin the distribution is flat-topped with re¬ 
spect to a Gaussian (with the effect becoming systemati¬ 
cally stronger when approaching the virial radius). This be¬ 
haviour is reflected in the decreasing kurtosis profile in the 
lower right panel of Fig. 0 The distribution shown in the 
second panel turns out to be quite well fitted by Tsallis dis¬ 
tribution f{v) tx [1 — (1 — with entropic 

index q = 0.65 and vo = 1178 km s“^ which we show with 
a solid line (see Lavagno et al. 1998; Hansen et al. 2004). 

The third panel shows the velocity distribution in the 
infall region: the mean velocity of the particles is now neg¬ 
ative, but the shape of the distribution is also altered. The 
distribution is more peaked than a Gaussian signifying a 
dominance of ordered flow over random motion of particles. 
The last panel illustrates the velocity distribution in the 
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Vj. [lO^m/s] 




Figure 4. The radial velocity distribution of the composite clus¬ 
ter at different distances. The histograms were made for radial 
bins of equal number of particles. The range of particle positions 
in each bin is written in the corresponding panel. On top of the 
histograms we plot with dashed lines the normalized Gaussian 
distributions with dispersions measured in a given bin. The solid 
line in the second panel shows a best-fitting Tsallis distribution 
with q = 0.65. 

outflow region, i.e. at distances larger than the turn-around 
radius. Here the distribution is also more peaked than the 
Gaussian but the mean velocity of particles is positive. 

We have therefore shown that, contrary to the still com¬ 
monly held belief, the velocity distribution in dark matter 
haloes is not Gaussian. Except for the very centre, it re¬ 
mains flat-topped inside the virial radius. The transition 
from the virialized to the non-virialized region is marked 
by a change in the velocity distribution from flat-topped 
to strongly peaked with respect to the Gaussian. We have 
demonstrated that the velocity moments in the virialized re¬ 


gion are quite well reproduced by the solutions of the Jeans 
equations. 

The Gaussian shape of the velocity distribution in the 
centre is expected since this kind of distribution is charac¬ 
teristic of structures with isothermal density profiles (r“^) 
and isotropic orbits (Binney & Tremaine 1987; Hansen et 
al. 2004). Hansen et al. also show that for isotropic orbits 
steepening of the density profile results in flattening of the 
velocity distribution, as we observe in our haloes despite the 
fact that their orbits are mildly radial. It would be very in¬ 
teresting to fully understand the relation between the shape 
of the distribution and the density profile and anisotropy of 
particle orbits. 
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